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Abstract 

We investigate the formation process of nonlinear vibrational modes rep- 
resenting broad H-bridge multi-site breathers in a DNA-shaped double 
strand. Within a network model of the double helix we take individual 
motions of the bases within the base pair plane into account. The result- 
ing H-bridge deformations may be asymmetric with respect to the helix 
axis. Furthermore the covalent bonds may be deformed distinctly in the 
two backbone strands. Unlike other authors that add different extra terms 
we limit the interaction to the hydrogen bonds within each base pair and 
the covalent bonds along each strand. In this way we intend to make 
apparent the effect of the characteristic helicoidal structure of DNA. We 
study the energy exchange processes related with the relaxation dynamics 
from a non- equilibrium conformation. It is demonstrated that the twist- 
opening relaxation dynamics of a radially distorted double helix attains 
an equilibrium regime characterized by a multi-site H-bridge breather. 

1 Introduction 

Studies of chemical and physical properties of DNA have attracted considerable 
interest among physicists and biologists because of their relevance for a variety of 
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biological processes, such as DNA transcription, gene expression and regulation 
and DNA replication 1 . In the context of the transcription process, the coding 
sequence on a DNA strand has to be made accessible to the RNA polymerase 
which necessitates that hydrogen bonds connecting the two strands have to be 
(temporarily) broken so that the DNA strands separate. This melted segment of 
DNA encloses 15 — 20 opened base pairs and is called the transcription bubble. 

The opening of the DNA molecule is a very complex process and a com- 
prehensive explanation of the actual mechanisms underlying the observed dy- 
namical processes is still to come. To tackle the problem of DNA dynamics 
nonlinear models of the DNA double helix have been proposed during the past 
two decades EJ'El- Different from molecular dynamics simulation methods, for 
which a great deal on the detailed reproduction of the molecular DNA structure 
is spent, the nonlinear models focus on the most relevant structural features only. 
The strength of the rather abstract nonlinear models of DNA lies then in the 
straightforward manner in which their soliton and breather solutions reproduce 
the prevalent dynamical behavior related with the strong energy localization 
and stable energy transport observed experimentally during the transcription 
process. 

To model the base pair opening in DNA, Peyrard and Bishop (PB) pro- 
posed a planar ladder-like model of DNA assigning each base pair a vertical 
inter-strand vibrational degree of freedom simulating the stretchings and com- 
pressions of the corresponding H-bridges [2] • The binding forces of the hydrogen 
bridges are described by a Morse potential. The bases itself are treated as point 
masses. Horizontally, the bases on the same strand are coupled via the cova- 
lent interaction described by harmonic potentials. The PB has been extensively 
studied and localized oscillating solutions (breathers) and moving localized ex- 
citations ^21 have been found reflecting successfully some typical features of the 
DNA opening dynamics such as the magnitude of the amplitudes and the time 
scale of the breathing of the 'bubble' occurring prior to thermal denaturation 
|2 ■ In addition, studies in the context of thermal denaturation have been based 
on the PB model incorporating a heat bath |14 | -|16 j . 

Typical for many biomolecules is their interplay between functional processes 
and structural transitions [Q. In this respect the transcription process relying 
on the opening of DNA represents a prominent example. Indeed, it has been 
found that the bubble formation is strongly correlated with twist deformations 
and local openings are always connected with a local untwist of the double helix 
[H] ■ [H] - - -"- n order to account for the helicoidal structure of DNA the PB model 
has been significantly extended by Barbi, Cocco and Peyrard (BCP) [Sj- In their 
model of a more realistic description of DNA, two degrees of freedom per base 
pair are introduced. There is a radial variable measuring the distance between 
two H-bridged bases along a line that connects them in the base plane being 
perpendicular to the helix axis. Further, the twist angle between this connecting 
line and a reference direction determines the orientation of the H-bridge. The 
twist-opening dynamics is described by radial breather solutions combined with 
kink-like solutions in the angular variable which have been constructed with the 
help of a multiple scale expansion technique 
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In the current study of the nonlinear dynamics related with the opening 
process in DNA our aim is twofold. First of all, we study the energy exchange 
processes and the relaxation dynamics in DNA molecules after their excitation 
into a non-equilibrium conformation. Such an investigation is associated with 
recent mechanical experiments performed with single DNA molecules |17 ] -|21 | 
forced away from their equilibrium conformations. After force applications en- 
ergy redistribution within the DNA molecule takes place such that a new equi- 
librium conformation is attained |22| . Furthermore, for strong enough radial 
forces the mechanical unzipping of single DNA molecules is achievable |2H|- We 
remark that for the current study our assumption is that the DNA molecules 
are imposed to not too strong external radial forces excluding the (immediate) 
separation of the strands. Secondly besides the relaxation dynamics within B- 
DNA whose two strands are supposed to have been pulled apart in a certain 
region, we direct our attention on the formation of oscillating bubbles occurring 
previous to denaturation. In particular, within a nonlinear model approach we 
aim for the creation of multi-site radial breathers, which, compared with the 
one-site breather solutions obtained in [HI, Eli provide stronger opening of an 
fairly extended segment of the DNA reproducing the experimentally observed 
oscillating bubbles involving 15 — 20 base pairs. 

The structure of the double helix of B-DNA is modeled by a steric network 
of oscillators in the frame of the base pair picture [2], [HI taking into account 
deformations of the hydrogen bond within a base pair and twist motions between 
adjacent base pairs. In augmentation of the oscillator model for the helicoidal 
DNA structure introduced in [Hj we allow for individual motions of the bases 
within the base pair plane such that the vibrations of the H-bridges are no longer 
exclusively symmetric with respect to the helix axis unlike the symmetric radial 
motions considered in [SJ. 

It is worth mentioning that when constructing helical variants of the Peyrard- 
Bishop model, the description of the hydrogen-bridge and covalent bonds is 
fairly clear: an asymmetric soft potential for the hydrogen bonds and harmonic 
terms for the covalent bond, in the hypothesis that the latter are more rigid 
and therefore the oscillations are small. It is not so evident how to model the 
stacking interaction in a straightforward manner. The physical origin of it lies 
in the bonds between overlapping bases, but these do not appear explicitly in 
the model. Different heuristic solutions have been tried as it will be commented 
in more detail later with interesting results. Here we focus on other aspects 
of the problem and take the helical shape of the DNA as given and add no 
extra term to maintain it. In this way we pretend to explore the effect of the 
helical shape in itself. The results are remarkable, the helical shape provides 
a change of the type of the effective coupling bringing about the possibility of 
stable multi-site H-bridge breathers, which have to be differentiated from broad 
one-site breathers. In this way much larger openings of the double strand are 
possible. 

The paper is organized as follows: In the second section we describe our 
extended steric network model for the helical structure of the double helix of 
B-DNA. The third section deals with the relaxation dynamics within B-DNA 
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molecules forced into locally distorted configurations and the fourth section is 
devoted to the construction of multi-site radial breathers. In the fifth section 
we discuss the linear modes of the breathers and their relation to the spectral 
features of the localized solutions. Finally, we present our conclusions. 

2 Oscillator model for the helical structure of 
DNA 

In our DNA model we incorporate the basic geometrical features of the DNA 
double helix structure which are essential to model the nonlinear dynamics of 
the twist-opening process. Similar to the approach in we treat the right- 
handed helical form of B-DNA, which is a polymeric molecule composed of two 
coiled strands of nucleotides forming a double helix, as a helical double-stranded 
oscillator system. The constituents of the latter represent the nucleotides which 
are regarded as single nondeformable entities. Thus, no inner dynamical degrees 
of freedom of the nucleotides are taken into account which is justified by the 
time scale separation between the small-amplitude and fast vibrational motions 
of the individual atoms and the slower and relatively large-amplitude motions 
of the atom groups constituting the nucleotides pQ . Concerning the structural 
components, each nucleotide is composed of a sugar, a phosphate and a base. 
The sugar-phosphate groups of neighboring nucleotides on the same strand are 
linked via covalent bonds establishing the interaction related with the rigid 
backbone to the strand. There is a base attached to every sugar. Since, for 
simplicity, we do not distinguish between the four different types of bases, the 
nucleotides are considered as identical objects of fixed mass. Two bases on 
opposite strands are linked via hydrogen bonds holding the two strands of DNA 
together. 

In the BCP twist-opening model, based on the base pair picture of DNA, 
the helicoidal structure of DNA has been conveniently described in a cylindrical 
reference system where each base pair possesses two degrees of freedom, namely 
a radial variable measuring the distance between a base and the central helix 
axis (viz., deformations of the corresponding H-bridge) and the angle with a 
reference axis in a plane perpendicular to the central axis which defines the 
twist of the helix 

Our model approach of the structural dynamics of the helicoidal DNA mod- 
els is inspired by the ones used in Ref. |5|-|10|. but differs from them in two 
important points: a) We release the constraint that the two bases in a base pair 
perform solely vibrational variations of the hydrogen bond length symmetric 
to the central axis along a line which connects the two bases crossing the cen- 
tral helix axis; b) Those models introduce for convenience extra potential terms 
which we discard because of the following reasons: On the one hand, not only 
they are sufficiently justified on physical grounds. On the other hand, we are 
interested in the effect of the helical structure on breather formation and try to 
keep the model as simple as possible to avoid masking the dynamical behavior. 
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These differences will bring about remarkable consequences on multibreathcr 
stability and on the phonon spectrum as will be shown in the next sections. 
We also neglect vertical movements as Ref. shows that they are not signifi- 
cant. Therefore, we extend the BCP model in the sense that we treat the two 
H-bridged strands individually and take individual two-dimensional equilibrium 
positions as well as displacement coordinates of a base within its base pair plane 
into account. As a consequence, the line connecting the two bases of a base pair, 
which is supposed to coincide with the orientation of the hydrogen bridge, does 
not necessarily intersects the central axis as distinct from the BCP model. In 
other words, for the equilibrium configuration of the double helix the two bases 
of a base pair may be positioned asymmetrically in radial direction and may be 
rotated also by different angles around the helix axis allowing for the simulation 
of polymorphic helical DNA matrices deviating from the perfectly regular equi- 
librium helix structure. In fact, real DNA molecules exhibit random structural 
imperfections of their equilibrium double helix caused by external and inter- 
nal influencing factors such as the distorting impact of the solvent environs. 
Moreover, with our approach the DNA lattice dynamics can be initialized with 
arbitrary excitation patterns connected with individual displacements of each 
base different from the symmetric H-bridge deformations discussed previously 
in 0. 

We describe the double helix structure in a Cartesian coordinate system 
whose z— axis corresponds to the central helix axis as sketched in Fig. ^ The 
base pairs are situated in planes perpendicular to the central helix axis and 
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the vertical distance between two consecutive planes is given by h. For the 
equilibrium configuration of B-DNA each base, with equilibrium coordinates 
\ and y£\, is rotated around the central axis by an ang le 6 { °\. The index 
pair [n, i) labels the n-th base on the i— th strand with i = 1, 2 and 1 ^ n ^ N , 
where AT is the number of base pairs considered. The equilibrium distance 
between two bases within a base pair, do, is determined by 

d, = \J{dl? + {dl)\ (1) 

where d n = x^\ — x^\ and d n = y n °\ — y n °\ are the projections of the line 
connecting the two bases on the x, y— axes of the coordinate system. The devi- 
ations d n from do through displacements, x n ,i and y n ,i, of the bases from their 
equilibrium positions, x^ \ and y n °\, are expressed as 

d n = \f (d n + x n ,i - x n . 2 f + ( d v n + y~i - y n . 2 f - do . (2) 
For later use we introduce the quantity 

O , d l + Vn.l ~ Vn,2 , „ 

n = arctan —y h 2 m ir , (3) 

dn "T 3*71,1 <Kn,2 

as the angle between the x— axis (as the reference direction) and the line connect- 
ing two (displaced) bases of a base pair measuring the alignment of the related 
H-bridgc, m being a integer to assure the monotonicity of 8 n with respect to n. 

The three-dimensional equilibrium distance between two adjacent bases on 
the same strand is given by 

io = ^ 9 +(*S-A,i) a + («23-A,i) a . w 

and deviations from l are determined by 



l n ,i = yh 2 + (L n l + X n a ~ X n -la) + ( L n i + y n>i - y n -i,i) - h , (5) 



model is then of the form 



-n,i ^n-i.i and L L = Wnj ~ fn-i,i- The Hamiltonian of our 



ff = E km + U h + U c , (6) 

where Uh and U c represent the potential energy part for the H-bond and the 
covalent bonds, respectively. The kinetic energy is determined by 



N 

m 



^-i EE [(*')' +(««) ! 



2m 

l= l,2 ra=l 

where m is the mass of a base and p^ 2 '' 1 denote the (a;, y)— component of the 
momentum of a base. 
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The vibronic potential Uh represents the deformation energy of the hydrogen 
bonds linking the two bases of the pair. The dynamical deviations from the 
equilibrium bond length d n (t) are supposed to evolve in a Morse potential and 
the potential energy related to the displaced hydrogen bonds is given by 

N 

Uh = D [exp(-ad„) -if , (8) 

n=l 

a^ 1 and D are respectively the width and the depth of the Morse potential 
well, corresponding the latter to the dissociation energy of the H-bond. In our 
simplified model of the DNA double helix we do not distinguish between the 
two different pairings in DNA, namely the G-C and the A-T pairs. The former 
pair involves three hydrogen bonds while the latter involves only two. 

Compared to the weak and flexible hydrogen bonds (with bond energies of 
the order of 0.04 — 0.3 eV) the covalent bonds between the sugar-phosphate 
groups of neighboring nucleotides on a strand are rather strong and rigid (with 
bond energies of the order of 2 — 10 eV). Therefore it is appropriate to treat 
the potential of the covalent bonds, simulated by elastic rods, in the harmonic 
approximation given by 

JV 

Po = ifEEC ( 9 ) 

i=l,2 n=l 

where K is the elasticity coefficient. 

With regard to other dynamical degrees of freedom we remark that longi- 
tudinal acoustic motions along the strands are significantly restrained by the 
rigidity of the phosphate backbone • Therefore we discard displacements of 
the bases in z— direction. Hence, the structural dynamics of base motions is 
restricted to the base planes. 

We will use the values of the parameters for B-DNA molecules given by 
(2,11,0: a = 4.45 A" 1 , r « 10 A, h = 3.4 A, Q = 36°, D = 0.04eV, K = 
1.0 eV A" 2 and M = 300 amu = 4.982 X 10~ 25 k g. 

With a suitable time scaling t — » Do? I rat one passes to a dimensionless 
formulation with quantities: 

~ _ ~(x) "n,i ~(y) "n,i fir\\ 

%n,i — M n,i i Un,i — OiUn,i > P n j , t P n i / „ J 

V mD yrnD 

K 

K = ' d n =ad n , r = ar , h = ah. (11) 

a z D 

In the following, we omit the tildes. 

The equations of motion are derived from the Hamiltonian (jSJ and read as 



x 



1 1 . i 



- J*) 



pZI . (12) 



Pn} = 2 [exp(-d„) - 1 ] exp(-d„) 
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2K 



, dl n .i dl n +i y 



9x T , 



iv) 



rn.i i 



P 



with the derivatives 
dd n 



2 [exp(-d n ) - 1] exp(-d„) 



dy n s 



2K 



&ln,i 



+ 1 



n+1,2 " 



dl n +i, 
dy n ,i 



dx n . 

dl n ,, 

3x n . 



L x -\- x n ^i — x n —\ i 

\J h2 + ( L n,i+ X n,i - X n-l,i) 2 + ( L' nl + y nA - y n _ i^) ' 



(13) 
(14) 

(15) 

(16) 
(17) 



and the equivalent expressions for dd n /dy n ,i and dl n ,i/dy n .i- 
It is useful to compare our model with previous ones: 

1. Refs. [HIE] add a three body curvature term Go {4> n +i + <t>n-i ~ 2 4> n ) 2 
with a quite large constant G = 50eV/A _1 , the elastic constant being 
K = leV/A 2 . The reason to introduce this term is to guarantee the 
correct helical shape instead of an alternatively possible zig-zag one. 

2. Refs. discard the term above and add a stacking term E exp(— b(r n + 
r n+ \ — 2 i?o)) (r n — r„_i) 2 , a extra degree of freedom, the axial displace- 
ment h with an elastic energy K [h n — h] 2 instead of the covalent harmonic 
bond along the strand. The stacking term is proposed to describe the 
shear force that opposes sliding motion of one base over another. The 
exponential attenuation becomes important for large opening at tempera- 
tures close to dcnaturation because the decrease of the molecular packing. 
The elastic energy with h < l guarantees the helical shape. The value 
of E — 4eV/A 2 is quite larger than K = 0.014 eV/A 2 , bringing about 
probably the dominance of the stacking term over covalent terms. 

3. Ref. jll| maintains the covalent energy term, discards the curvature term 
and adds a stacking term but without the attenuation factor, as he does 
not consider high temperatures. (See the discussion of the different terms 
in this reference). 

4. Ref. |10j uses the harmonic stacking term on the radial variables, the 
curvature term, the covalent bond term and the vertical elastic term. 

In our work we do not introduce all these terms, because we are interested 
in the effect of the helical shape in itself on breather formation, and they are 
not deduced from first physical principles. This makes it incomplete, as there 
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are, certainly, forces that maintain the helical shape, but we think that the 
introduction of extra energy terms would mask it. 

We think that an appropriate deduction of the form of the stacking energy 
terms is still to come. 

We use the following values of the scaled parameters K = 0.683, r$ — 44.50, 
h = 15.13 and lo — 31.39. One time unit of the scaled time corresponds to 
0.198ps of the physical time. 

3 Energy redistribution, relaxation dynamics and 
breather formation 

In the following we consider the energy exchange process between the radial 
and torsional degrees of freedom when a confined region of the DNA molecule 
gets deformed from its equilibrium configuration. Nowadays, there exist sophis- 
ticated experimental techniques for the selective excitation of DNA in single 
molecule experiments (see, e.g., and during the last years, mechanical 

properties of DNA molecules have received a lot of attention ^Tj [22- Several 
force measurements were performed on single DNA molecules to examine, e.g., 
the longitudinal extension ^B] and the twist elasticity of DNA molecules \17\ . 
Furthermore, the opening of two-stranded DNA molecules was mechanically 
forced by pulling apart the two strands of a DNA double helix [2U] , 

Besides the study of the relaxation process in deformed DNA molecules, our 
aim is also to create spatially extended H-bridgc breather solutions, reproducing 
the oscillating 'bubbles' observed for the DNA-opening process, which extend 
over 15 — 20 base pairs [S]. To this end we excited initially twenty consecutive 
lattice sites in the center of the DNA lattice assuming that the DNA molecule 
experienced deformations in radial direction. For the numerical simulation we 
elongated each of the associated twenty hydrogen bonds from its equilibrium 
length by displacing the Morse oscillators out of their rest positions accordingly. 
Provided these elongations are aligned solely along the equilibrium orientation of 
the hydrogen bonds no twist deformations occur (yet) . Naturally, the distortion 
of a hydrogen bond is connected with a deformation of the covalent bonds of 
the phosphate backbone in the neighborhood of the corresponding base pair. 
Nevertheless, for the starting excitation pattern, with radial amplitudes in the 
range of rf{„ c } = (0.1 — 0.35) A, one finds only small deformations of the covalent 
bonds. The set {n c } labels the indices of the twenty excited lattice sites in the 
central region. Compared to the equilibrium DNA conformation the deformed 
one possesses an amount of potential energy increased by the deformation energy 
(also referred to as excitation energy). Furthermore, the energy contained in 
the elongated hydrogen bonds, Uh, is exceedingly larger than those contained 
in the deformed covalent bonds, U c , and the ratio is typically of the order of 
U c /U h < 0.15. 

We integrated the set of coupled equations (|T^|) - ((r3|l with a fourth-order 
Runge-Kutta method while the accuracy of the computation was checked through 
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monitoring the conservation of the total energy. For the simulation the DNA 
lattice consists of 400 sites and open boundary conditions were imposed. (Wc 
remark that, provided the system size is sufficiently large, our findings are in- 
sensitive to further enlargement of the system.) First of all, we report on the 
results obtained for the regular equilibrium configuration of the DNA molecule 
for which the bases of a pair are rotated around the axis by the same angle 
6 n ,i = On.2 +7T = Or? = ti,0q. Consequently, the equilibrium twist angle between 
two consecutive base pairs is 9o and all bases possess an equal radial distance 
do/2 from the central axis. In Fig. 0-a we depict the spatio-temporal evolu- 
tion of the distance, d n (t), between two bases of a base pair, which measures 
the variation of the length of the corresponding hydrogen bond. Initially, the 
excitation energy is divided evenly between the twenty lattice sites so that the 
initial distance profile, d n (0), is rectangularly shaped. In the current regular 
case, the two bases of a base pair move symmetrically with respect to the cen- 
tral axis in radial direction due to the asymmetric choice of initial conditions, 
(^n,i(0), 2/n,i(0)) = — (x„ i 2(0), y n ,2(0)). Thus, the associated displacement coor- 
dinates perform out-of-phase oscillations (x n ^i(t), y n ,i(t)) = —{x n ,2{t),y n ,2{t))- 
In addition, the two bases of a pair possess symmetric angular dynamical dis- 
placements, 6 rh \(t) = 8 n .2(t) + ir = 9 n (t), and the line connecting the two 
(elongated) bases in the base plane always intersects the central axis. The dis- 
tance variable d n (t) is hereafter also referred to as the radial variable because 
d n (t)/2 represents actually the local helix radius. The corresponding angular 
displacement pattern, 9 n (t) — 9n \ is shown in Fig.|2-b. 

As becomes evident from Fig.|2-a the vast majority of the excitation energy 
remains in the initially excited central region and for the radial motion a multi- 
site breather develops for which the stretching of a base pair distance is larger 
than the compression characteristic for the evolution in a Morse potential (see 
also jHJ). Moreover, the resulting spatially extended radial breather involves 
all of the initially excited twenty oscillators evolving with practically equal am- 
plitudes. To either side of the central region the amplitude pattern abruptly 
goes to zero giving the multi-site breather an almost rectangular spatial pro- 
file nearly matching the initial profile. The maintenance of the broad localized 
radial shape has to be distinguished from the dynamical behavior reported in 
[7] where, in the context of the BCP model, for similar initial conditions in the 
radial stretchings in a fairly broad region a merging of the radial components 
is observed in the course of time which results eventually in a narrow radial 
pattern localized at a single site only. 

For the associated angle deformations 9 n (t) — 0n \ emerging from overall 
initial zero amplitudes, an asymmetric pattern generates rather rapidly in the 
initially excited region. It holds that the angles from the left (right) of the two 
central lattice sites (base pairs) decrease (increase) steadily and, the further 
apart a base pair is from the central ones, the stronger its angular displace- 
ment. Eventually, a local kink-like pattern in the angular lattice 9 n {t) — 6n 
is created. The horizontal plateaus of the kink-pattern extend continuously in 
either direction away from the central base pairs so that in the course of time 
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Figure 2: Relaxation dynamics: For a proper illustration of the extension of the 
radial breather and the forming kink-like angular structure an inner segment 
of the DNA lattice consisting of 200 sites is shown. Initial conditions: rf{„ c j = 
0.1 12 A. (a) Radial breather, (b) Angular deformation pattern. Development of 
the kink-like structure. 
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more and more base pairs become subject to angle deformations leading to a 
progressive untwisting of the central part of the helix j23| ■ This untwisting of 
the helicoidal helix structure results from the coupling between the radial and 
the torsional degrees of freedom due to geometrical constraints and is typical 
for the DNA opening dynamics UJ. In contrast to the periodically oscillating 
pattern of the radial variable d n (t), corresponding to alternate stretchings and 
compressions of the hydrogen bonds, the torsional deformations 6 n (t) — 9n^ ad- 
just to a static deformation pattern. For a direct comparison of the extension of 
the radial to the angular helix deformations, the angular component, expressed 
in radians, has to be multiplied by r « 10 A to estimate the associated length 
scale. We conclude that the degree of the deformations in the angular direction 
are comparable with those in the radial direction. This has to be distinguished 
from the observations made for the evolution of narrow radial breathers for 
which the dynamics of the opening process is dominated by the displacements 
of the bases in radial direction whose amplitudes are typically larger than those 
of the angular displacements by at least one order of magnitude 

For further illustration of the energy sharing phenomena accompanying the 
breather formation process, we monitored the time-evolution of the energy con- 
tained in the initially excited 20 base pairs, which is determined by 



Ecenterij^) — ^ ^ 
{™c} 



+ U h (d n )\. (18) 



The temporal behavior of E center {t) is represented by the graph with label (1) 
in Fig. [fj] We observe that in an early phase, lasting for approximately forty 
time units (~ 8ps), a small energy loss takes place and around 1% of the 
excitation energy flows from the initially excited central region into wider parts 
of the remainder of the DNA lattice. Afterwards the amplitude of the energy 
E ce nter(t) fulfills slow oscillations around a mean value. (That this mean value 
of the energy E center (t) is weakly descending as time progresses indicates further 
small energy emission from the central region into the rest of the DNA lattice.) 
The energy emitted from the central region is then mainly redistributed into 
the angular deformation kink that travels through the DNA lattice in agreement 
with the formation of the angular pattern represented in Fig.|2|-b. Nonetheless, a 
tiny part of the emission energy is radiated also into the hydrogen bonds outside 
the central region where it is converted in radial phonons moving uniformly 
towards the ends of the lattice as just about visible in Fig. EJ-a. 

The Fig.|31shows the distribution of the energy, between the hydrogen bonds 
and the covalent bonds, which has been deposited from the central region into 
the rest of the DNA lattice. Plotted are the time-evolution of the potential 
energy contained in the displaced hydrogen bonds 

U c h (t) = £ (eM-drM) ~ I) 2 , (19) 
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Figure 3: Relaxation dynamics: The distribution of potential energy between 
the hydrogen bonds and the covalent bonds outside the central region of the 
DNA lattice. 



and the potential energy content of the covalent bonds 

u c c {t) = K 

n, i=l,2 



(20) 



respectively, where n <G [1, N]\{n c }. After the initial phase of energy absorption 
a quasi-equilibrium regime is reached and the partial potential energies oscillate 
around mean values of 0.0007 eV and 0.0003 eV for the deformations of the 
covalent bonds and the hydrogen bonds, respectively. 

Complementary, the energy sharing between the deformed covalent and hy- 
drogen bonds within the initially excited region is illustrated in Fig. 0] During 
a transient phase of internal energy redistribution the average of the covalent 
bond energy grows slightly on expense of the hydrogen bond energy. However, 
the energy migration is not very pronounced. Eventually, the time evolution of 
each of the two partial potential energies is characterized by periodic oscillations 
around a (constant) mean value reflecting the attainment of a steady equilib- 
rium regime in the central lattice region. Note that the initial drastic difference 
in the energy contained in the stretched hydrogen bonds, Uh{0) = 0.124 eV, and 
the covalent bonds, U c (0) = 0.018 eV, is retained throughout the time-evolution 
of the respective average partial potential energy. 

In conclusion, we found that out of an initial non-equilibrium situation, for 
which the hydrogen bonds in a fairly broad but confined region of the cen- 
tral part of the DNA lattice have been stretched away from their rest lengths, 
(weak) energy redistribution processes set in, such that small amounts of energy 
migrate into the remainder of the DNA lattice. During the relaxation process 
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Figure 4: Relaxation dynamics: The energy sharing between the deformed cova- 
lent bonds (upper graph) and hydrogen bonds (lower graph) within the initially 
excited region. 

towards an equilibrium of energy balance between the radial and torsional com- 
ponents the angular displacement variables adopt a static kink-like structure. 
Correspondingly, due to the torsional deformations induced by the opening of 
the base pairs a local unwinding of the helix develops. Significantly, almost 
the entire initial excitation energy stays localized in the hydrogen bonds of the 
central region and a multi-site breather with nearly rectangular profile develops 
in the radial displacement variables. 

Besides the initial pattern of broad radial inflation we imposed initially also 
radial compressions in a fairly broad region of the DNA lattice giving also rise to 
multi-site breathers corresponding to alternating compressions and stretchings 
of the H-bonds. 

Finally, regarding asymmetric initial conditions we mimicked e.g. the im- 
pact of a local 'pushing' force in cause of which some bases are brought closer 
to the helix axis. In Fig. [3] we display the spatio-temporal evolution of the ra- 
dial deformations of the H-bonds when initially three nucleotides at only one 
strand were displaced from their equilibrium positions. Apparently, even for 
such asymmetric situations H-bond breathers result characterized by deforma- 
tions of the hydrogen bonds which are no longer symmetric with respect to the 
central axis. The conclusion is that the multibreathers are not destroyed by 
interaction with the modes of the centers of mass . Other numbers of excited 
base pairs bring about similar results. 
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Figure 5: H-bond breather formation for asymmetric initial conditions. Three 
consecutive sites on one strand possess initially diminished distance to the helix 
axis. 

4 The relaxation concept and breather construc- 
tion 

We exploit now the findings of the preceding section for the construction of 
multi-site breathers in DNA. For this goal we initialize the DNA lattice as 
described above, that is twenty Morse oscillators in the central part of the 
DNA lattice are elongated from their equilibrium lengths. Subsequently we let 
the DNA lattice dynamics relax towards an equilibrium regime during a time 
interval of 100 time units. We recall that for such times the energy redistribution 
between the hydrogen bonds and the covalent bonds and the central region 
and the remainder of the DNA lattice, respectively has proceeded so far, that 
virtually negligible alterations of the partial energy contents occur for later 
times. Notably, a spatially extended radial breather forms on such a time scale. 
Furthermore, due the geometrical constraint of the helicoidal structure a local 
untwisting of the double helix arises. The corresponding kink-like pattern in 
the angle variables encompasses almost the whole lattice, apart from the small 
parts of the lattice near the boundaries which have not yet experienced angle 
deformations (cf. Fig. |5J-b). 

We start now a iteration scheme, for which the quasi-equilibrium solution 
resulting at the end of the simulation time, viz. the displacements (x ni i(100), 
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J/n,i(100)), is taken as the new initial condition for the re-iterated DNA lattice 
dynamics. In order to attain an equilibrium regime of the DNA lattice dy- 
namics the moving radial phonons, carrying the excess energy emitted from the 
central region, have to be eliminated. In the primary iteration step the radial 
phonons transport energy merely of the order of 10 -4 eV. With each further 
iteration cycle the amount of radial phonon energy, still ejected from the central 
lattice region, further diminishes. Nevertheless, to accomplish faster relaxation 
we impose absorbing boundary conditions to the DNA lattice so that the radial 
phonons are removed, once they arrive at the ends of the lattice. In this iterative 
manner the energy redistribution from the central region to the remainder of the 
DNA lattice is more and more reduced and the dynamics relaxes towards a solu- 
tion of improved energy balance compared to the previous step. To illustrate the 
latter fact, we display in Fig. the temporal evolution of the energy, E center (t), 
contained in the central part of the lattice which bears the radial breather. Ap- 
parently, convergence of the iteration procedure is achieved already after three 
iteration steps. In Fig. 0we depict the resulting radial breather solution. This 
breather oscillates with a period of = 0.7ps which is a realistic value for the 
time scale of experimentally observed vibrational modes occurring prior to the 
opening process in DNA melting More importantly, the broad localized 

radial excitation pattern comprises twenty base pairs and the maximal elonga- 
tion of a base pair radius from rest length is about 0.122 A . Thus, we are able 
to create long-lived spatially extended radial breathers of fairly large amplitude 
in DNA representing the oscillating 'bubbles' as the precursor to thermal de- 
naturation. Finally, we remark that, depending on the choice of the width of 
the initial excitation pattern, with our relaxation-method radial DNA breathers 
of varying extension can be generated. These breathers range up from nar- 
row ones (effectively one-site) up to very broad breathers involving up to even 
fifty lattice sites. As the choice of the initial excitation pattern is further con- 
cerned, we tried, besides the radial rectangular profiles, also other profiles such 
as bell-shaped ones given by d{„ c }(0) = A / cosh[/3 (n — n )]. A , (3 and n are 
constants fixing the amplitude, the spatial extension and the center of the start- 
ing profile, respectively. In the outcome of the relaxation dynamics we observed 
bell-shaped radial breathers nearly preserving the initial profile analogous to the 
results reported earlier in this manuscript for the broad breathers. The same 
holds when additionally to the radial displacements also initial kink-like twist 
deformations are considered. 

For a more elaborate study of nonlinear vibrational modes in DNA the ef- 
fects of external influences have to be included. Real DNA molecules exhibit 
random structural imperfections of their double helix caused, e.g., by the de- 
forming impact of the chemical surroundings when DNA gets buffeted by water 
molecules. Another source for structural irregularities originates from the ran- 
dom base sequence of the genetic code. In addition, the varying hydrophobic 
potential of the base pair interactions depending on the ambient aqueous solvent 
may leave the helix structure in an irregularly distorted shape. Accordingly, for 
an improvement of our model, irregularity effects can be mimicked by taking 
into account structural disorder so that the resulting irregular helical DNA ma- 
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Figure 6: Breather construction: The temporal evolution of the energy con- 
tained in the central part of the lattice for the first, the second and the third 
step of the iteration procedure, respectively, as indicated in the plot. The ener- 
gies are normalized in units of the initial energy of the first iteration E center (0). 



trix deviates from the perfectly regular helix structure. (The formerly discussed 
ordered, regular structure arises, for example, for synthetically produced DNA 
molecules consisting of a single type of base pairs, for instance, poly(G)-poly(C) 
DNA polymers, surrounded by vacuum.) 

To be precise, we consider randomly distributed equilibrium coordinates, 
xfj and y^J, of the bases with i^ ) - x { °\ e [-Ax n>i , Ax nii ] and y^J - y^J € 
[— Ay„.i, Ay Ui i\. The intervals of the deviations Ax and Ay range up to 5% of 

the respective rest value x^\ and y£\ of the corresponding regular structure. 

In Fig. [S] we show a representative case for the DNA lattice dynamics for 
one realization of structural disorder. In addition to the disordered helicoidal 
structure also random amplitudes of the initial values for the distance displace- 
ments, d{ n }(0), in the central region have been chosen. Apparently, the H-bridge 
breather is robust sustaining the influence of disorder and virtually resembles 
the behavior of the ordered case (compare Fig.|5|-a). 

On the other hand, when the amount of excitation energy placed in the cen- 
tral lattice sites, exceeds a critical value, E Uc > 0.6 eV , or equivalently when the 
hydrogen bonds get elongated too far, we do not longer observe the relaxation 
towards a stable multi-site radial breather as in the previous (lower-energy) 
cases. Instead directed flow of energy into one of the initially excited lattice 
sites takes place related with the simultaneous energy depletion of the remain- 
ing lattice sites. In this manner, the transferred energy is accumulated in the 
corresponding hydrogen bond of the energy-gaining site with the result that the 
radial elongation continually grows. Eventually, the corresponding hydrogen 
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Figure 7: Breather construction: The final radial breather solution obtained 
after three iteration steps of the relaxation concept. 
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Figure 8: Breather construction: The persisting radial breather for a disordered 
helicoidal equilibrium structure (see text). The initial values for the irregular 
distance displacements in the central region, d{ nc }, are randomly distributed 
according to d{„ c } — di n \ G [—Ac?, Ad]. The interval of the deviations Ad 
ranges up to 5% of the respective value d{„ c } of the regular case depicted in 
Fig.0 

bridge is stretched so far that it breaks up. This energy concentration onto a 
single site and the succeeding bond breaking process is illustrated in Fig. [5] 

Conversely, even for as tiny radial displacement amplitudes as d{ n y(0) ~ 
10~ 5 A, corresponding to a lattice excitation energy of the order of 1CP 7 eV, the 
energy stays localized at the twenty initially lattice sites and a radial multi-site 
breather is formed. This demonstrates impressively that DNA possesses very 
efficient energy storing abilities. It should be mentioned, that usually, in studies 
of nonlinear oscillator networks one observes that for undercritical amplitudes, 
being tantamount to a low degree of inherent nonlincarity, an initially local- 
ized structure decays rather rapidly and the excitation energy gets uniformly 
distributed over the whole lattice |2fi| . 

It should be noted that Campa pHQ also shows the existence of broad bubbles 
that move along the molecule. The differences between his model and ours are 
commented in the previous section. The results are quite different. The Campa 
bubbles have no oscillation and a bell profile, while our multibreathers oscillate 
and have a quasi rectangular profile, involving a larger number of base pairs in a 
similar way. The oscillation is a major weakness of our bubbles as they can only 
be considered as precursors of transcription or denaturation. They have some 
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Figure 9: Concentration of the excitation energy on a single site of the DNA 
lattice and successive bond breaking. Initial conditions corresponding rf{ n< .} = 
0.75 AA corresponding to an amount of excitation energy of 0.646 eV. 



advantage, however, which are their quasi homogenous amplitude and that they 
are more generic, being the consequence of the helical shape, and not relying on 
evolving in the flat part of the Morse potential. 



5 Linear modes 

In this section we consider the linear modes of the system and their influence 
on the behavior and stability of the system. The linear equations corresponding 
to Eqs. (|12H15|1 for the homogeneous and symmetric system are given by: 

x n , i = (-l) i 2(<%Ax n +d»Ay)?Z (21) 

a 



2K&[L* .Sx nii + LlJy^] - ^hl[L^Sx n+1 , + Ll +1 Jy n+1 ,] ) , 

d y 

V n j = (-iy2(<%Ax n + d»*y)-£, (22 

"o 

L y L y 
-2K(-£±[L% ti 5x n ,i + L v ni 8y n ,i) r ^^[L x n+1 Jx n+1 , l + L v n+1 Jy n+1: ,] ) 

where Ax n = x n<1 - x„ j2 , Ay„ = y n ^ - y n<2 , Sx n = x Uti - x n -i ti and Sy n = 
Un,i — Un—i,i- In the symmetric equilibrium state it holds that L^ L 2 = — l = 
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L o cos(0i o) ), L v n 2 = -L y n l = L sin(<W > ), dn,x = do cos(6>i 0) ), d n , y = d a sin6>l°\ 
with L = ((L*) 2 + (L^) 2 ) 1 / 2 = c?o sin(0 o /2), being the equilibrium horizontal 
distance between nucleotides along one strand; (jffl being the angular cylindrical 
coordinate of the vector that goes from nucleotide n — 1 to n along strand 1 at 
equilibrium, i.e., 4>n^ = On — Oq/2 + ir/2; and 9^ — jiOq, as defined in Eq.|3J 
at equilibrium. 

The natural coordinate system for each base pair is given by the projections 
on the unit vectors e r = (cos(0« , ), sin(6*i°'), 0) and eg = (— sin(0^), cos(0„ ), 0), 
and the natural variables are the relative displacements of the base pairs Ar„ = 
(Ax n , Ay n , 0) = (x n ,i — X n> 2,y n ,i ~ 2/n,2)0), and the positions of their centers 
of mass in the horizontal plane: r n — (x n , y n ,0) — ((i„.i + x n ,2)/2, {Un.i + 
2/n,2)/2, 0). Therefore, the natural coordinates are: 

r n = Ar„ • e r = cos(e^ ] ) Ax n + sin(6>W) Ay n , 
Ar n ■ eg 



a 



l(-sin(0f )Az„+cos(0f )Ay n ) , 



do d Q 

r n = r„ • e r = cos(6^)W n + sin(6»^ 0) )y n , 

s n = r n -eg = -sin(9 ( V)x n +cos(9^)y n . (23) 

To first order in the perturbations {£ n ,i},{yn,i}, T n ~ d n represent the stretch- 
ings of the base pairs; a n ~ arctan[(d^ + Ay n )/(d* + Ax n )] — 6*1°^ (mod 7r) , 
the relative twist angle with respect to the equilibrium position. r n and s n 
are the displacements of the base pairs centers of mass in the directions of the 
bonds or orthogonal to them respectively. In terms of these variables, the linear 
dynamical equations l|21H22|l can be written as: 

- • 2 1 \ /v, - , - - i ■ r®\ ^0w_ _ x 

r n = -fJ> sin {—) (2r„ + r n+1 + r n -i) - M sm(y) cos(y) {s n +i - s n -i) 

s n = (i sm(y)cos(y) (r„+i - r„_i) - /i cos (—J (2s„ - s n+ i - s n _ij , 

r„ = -w 2 r„ 

- /j sin (— ) (2r„ +r n+1 +r n _ 1 ) - fi sin(— -)cos(— Jdo (a n+ i - a n -i 

doan = M s in (y)cos(y) (r„+i - r„_i) - ,u cos (y)d (2a„ - - 



with /i = 2K(L /l ) 2 . 

The first two equations describe the motion of the helix as a whole, their 
corresponding linear modes frequencies will be referred to as the cm (center of 
mass) branch. The last two equations represent the dynamics of the base pairs 
and lead to two branches, referred to as the bp (base pair) acoustic and optical 
branches. Note that the cm equations are independent of the bp equations and 
so are their modes. 
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The frequencies of the cm branch are given by: 

lo 2 = 2/i [ 1 - cos{9 )cos(q)} , (25) 

and a branch of zeros, which represents the possible motions of the helix to 
other equilibrium positions. 
The bp branches are: 

<4 = 2 ( w ° + 2 M I 1 ~ cos ( 6, °) cos («)]) 

± + 2 [1 - cos(6» ) cos(q)]) 2 -8/jw 2 cos 2 (6» /2) sin 2 (g/2) .(26) 

The sign + and — correspond to the bp optical and acoustic branches respec- 
tively. 

Fig. HUI -a shows the dependance of the phonon spectrum on the coupling 
parameter K and Fig.llUhb depicts the three branches for the value of K = 0.683 
used in this article. It can be seen that there is a phonon gap for values of 
K < K c with K c ~ 0.7. The Fourier spectrum of the multi-breathers found 
consists of a single frequency tOb slightly below loq — 2, the frequency of the 
isolatedly vibrating base pairs, and its harmonics. Furthermore, the modes in 
the neighborhood of lo^ have wave vector q = tt and therefore they are not much 
excited by the perturbation of the radial variables. 

A remarkable fact of the phonon spectrum for K < K c is the inversion of 
the optical spectrum with respect to the planar PB model ,2,. In the latter, 
the only variables are the distances between bases within each base pair, and 
the coupling proposed is a standard attractive one, which has been made more 
complex afterwards |^|[2E|, but without consequences on the phonon spectrum. 
The dynamical equations are given by: 

u n + V'{u n ) + e(2u n - u n -i - u n+1 ) = , (27) 

being V the Morse potential with V"(0) = ui 2 , > and e > 0. The dispersion 
relation is given by lo 2 = lo 2 + 4e sin 2 (q/2). There are some consequences for 
this system [23 123 : a ) the bottom (linear) mode has q = 0, and the top 
one has q = tt; b) one-site breathers are stable; c) the tails of a breather 
or multibreather consists of in-phase oscillators; d) multibreathers with all the 
oscillators in phase (derived from q = 0) are unstable and with all the oscillators 
out of phase (derived from q — tt) are stable |^] ESI IS OS]. If the value of e 
were negative, the conclusions would be reversed. 

In the model proposed in this article considering the helical shape of the 
DNA molecule, the optical spectrum is inverted as can be seen in Eqs. [21 
Due to the degeneracies of the model we have not been able to obtain exact 
breathers, but the Fourier spectrum just shows a frequency below the optical 
branch and its harmonics, which means that they are standard breathers in 
the usual meaning, i.e., localized, nonlinear, periodic oscillations. The tails of 
the multibreathers are, indeed, of the q — n type, according to this inversion. 
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Figure 10: Spectral analysis of the linear system, (a) Dependence of the phonon 
spectrum on the longitudinal elastic constant K. For values of K below K c ~ 0.7 
there exist a gap below the optical branch where the frequencies of the multi- 
site breathers lie. (b) Frequencies of the three phonon branches with respect to 
their wave numbers q. The it modes are the ones on the vicinity of the multi-site 
breather frequency. 



However, the transcendental consequence arising from the inversion of the ef- 
fective type of coupling is that the multibreathers are stable. This is a result 
that should not be overlooked: the helical shape of the DNA molecule provides 
a means for stable, large amplitude oscillations of base pairs groups, i.e, for the 
formation of the denaturation bubble. 



6 Summary 

In the present work we have studied the formation process of broad H-bridge 
breathers in DNA molecules. The double helix formation of DNA has been 
described in terms of an oscillator model relying on the base-pair picture. In 
this context each base pair possesses four degrees of freedom, namely a radial, 
an angular one, and two motions for the base pairs centers of mass. For an 
extension of the BCP model we have treated the two strands individually 
allowing for asymmetric vibrations of the H-bridges. As for the simplifications 
of our DNA model, we remark that we have not distinguished between the four 



23 



different base types and therefore, have treated each base as a single entity of 
fixed mass. Apart from this, we neglected also the different types of hydrogen 
bonds between the two different pairings in DNA, namely the bridging of the 
G-C and the A-T pairs by three and two hydrogen bonds, respectively. 

We concentrate in this paper on the effect of the helical shape of the DNA 
molecule, avoiding the introduction of other potential terms to make it apparent. 
The result is that this shape brings about a inversion of the effective type of 
coupling, with a inversion of the phonon spectrum, that allows for the existence 
of stable multi-site breathers. 

We have demonstrated that in the course of the relaxation dynamics of DNA 
molecules, which have been displaced in a certain region from their equilibrium 
configuration, an equilibrium state is attained. More precisely, starting from 
a radially distorted configuration involving several consecutive base pairs the 
twist-opening dynamics relaxes towards a localized state. The latter is built up 
from a multi-site H-bridge breather in combination with a static kink-like profile 
of the angular variables connected with the untwist of the helix. Excess energy, 
not to be contained in the radial breather, is emitted from the initially excited 
region in the form of radial phonons approaching either ends of the DNA lattice. 
The relaxation process and the attainment of a stable breather regime has also 
been observed for asymmetric initial conditions when e.g. bases on only one of 
the two strands are displaced from their equilibrium positions. 

We have exploited then the dynamical approach of an equilibrium state for 
the construction of breather states of the DNA. To this aim we have established 
an iteration scheme, such that, with increasing steps of applied iterations the 
dynamics gets closer to an equilibrium regime supporting a radial breather to- 
gether with a kink-like pattern in the angular components. It should be stressed 
that the obtained broad H-bridge breathers reproduce realistically the exten- 
sion of the oscillating bubbles observed in DNA. However, the almost rectan- 
gular shape of our H-bridge breather solution has to be distinguished from the 
single-site radial breather being reminiscent of an envelope soliton solution with 
a half-width of twenty base pairs obtained in the context of the BCP model in 
Ref. .3 ■ The multi-site breather with all of its constituents performing in-phase 
oscillations of equal amplitude renders all of the involved base pairs accessible 
to the functional process (e.g. the transcription) on an equal footing in contrast 
to the uneven and less efficient separation of the two strands associated with 
the bell- shaped amplitude pattern when the breather is centered at a single site 
so that merely the half width of the exponentially localized profile comprises 
15 — 20 base pairs. 

Our iteration scheme for the construction of H-bridge breathers is applicable 
so long as the amplitudes of the radial distortions do not exceed a critical value. 
However, for overcritically large amplitude (or equivalcntly, too large amount 
of excitation energy) the dynamics does no longer relax onto a broad H-bridge 
breather. Remarkably, one rather observes the directed flow of excitation energy 
in a single H-bond with the consequence that this bond may even break up. 

Finally, we note that with view to the role of structural disorder, we have 
found that the broad H-bond breathers sustain moderate amount of randomness 
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in the arrangement of the equilibrium positions of the bases as well as in the 
initial local displacement patterns. 
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